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Abstract 

The critical boundaries separating ordered from chaotic behavior in randomly 
wired S-state networks are calculated. These networks are a natural generaliza- 
tion of random Boolean nets and are proposed as on extended approach to genetic 
regulatory systems, sets of cells in different states or collectives of agents engaged 
into a set of S possible tasks. A order parameter for the transition is computed 
and analysed. The relevance of these networks to biology, their relationships with 
standard cellular automata and possible extensions are outlined. 
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1 Introduction 



Cellular automata (CA) [1] and random Boolean networks (RBN) [2] are two key ap- 
proaches in our understanding of complexity [3,4]. In both cases, discrete-time and 
discrete-space are used, together with a finite number of states. For CA, we have a 
dynamical evolution defined (for the one-dimensional case with a (2r -|- l)-neighborhood) 
by means of: 

ai{t + l) = (7j_r(i),(7i_^+i(i), ...,ai+r{t) ] (1) 

where i — 1,...,N and ai{t) e E = {0,1,. ..,5" — 1}. Here $ is a given (previously 
defined) rule [1]. Identical rules are defined for each {2r + l)-neighborhood. In contrast 
with this homogenity in neighbors and rules, RBN's are a particularly interesting class of 
dynamical systems which shares some properties with CA models, but where randomness 
is introduced at several levels. Now the dynamics is given by: 

ai{t + l)^Ai [ai,{t),ai,{t),...,ai^{t)] (2) 

where i — 1,...,N and ai{t) e S = {0,1}. Each automaton is randomly connected 
with exactly K others wich send inputs to him. Here Aj is a Boolean function also 
randomly choosen from a set J-'k of all possible JC-inputs Boolean functions. In spite of 
the random selection both of neighbors and Boolean functions, it is well known that a 
critical connectivity Kc exists where "spontaneous order crystallizes" [2]. In Kc a small 
number of attractors (~ 0{\/N)) is observed which show high homeostatic stability (i.e. 
high stability against minimal perturbations in single elements or Boolean functions) , and 
low reachability among different attractors. These properties are consistent with some 
observations of the genome organization [2,5]. A considerable effort has been dedicated 
over the last years to the analysis (both theoretical and computational) of the scaling 
behavior of the numbers and size of cycles close to the critical boundary [5-7]. 

Alternative models to this problem are provided by neural network- like aproaches [8] . 
In this case we have a dynamical system given by the set: 

K-l 

ai{t + 1) = ^Juaiit) + ^ JiMi)'Jj,ii){t) + 9i] (3) 

z=i 

where <l>(x) is a sigmoidal function which can be assimptotically approched by a sign 
function, $(a;) = sgn{x) which is sgn{x) = if a; < and sgn{x) = 1 otherwise. 
Not surpringly, neural networks with random connectivities show some of the properties 
displayed by RBN. Two well-defined phases are also described, which also depend upon 
the input connectivity [9]. 

Though the picture of on-off genes is a simplified one, in some cases it can be justified 
from general arguments based on mutual inhibition among genes [10]. Let us consider 
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two given genes which interact among them by means of two given proteins, whose con- 
centrations are indicated by the pair {x, y) . This can be used as a toy model of some well 
known properties of the life cycle of viruses (like the A-Phage [3,10]). A theoretical model 
of mutual inhibition is given by the following couple of cross-inhibitory equations: 

^ = X U.a) 

^^^^ y {A.h) 

where x,y & R+. Here n e N is a measure of the strength of the interaction. For ji — 1/2, 
the stability of the fixed point P* — (1/2, 1/2) is given by the eigenvalues of the Jacobi 
matrix: 

L,{n^{zl If) (5) 

which leads to stability forn < 2 and to instability otherwise. For n> ric — ^, P* becomes 
a saddle point and two new stable attractors {P+,P_} are formed (stable nodes). For 
strong interactions, one of the concentrations dominate over the other and, basically, the 
two new attractors are simply P+ pa (1, 0) and P_ ~ (0, 1). So the result of this interaction 
between continuous quantities leads to a basically binary outcome [11]. 

But gene activity can lead to more complex combinations (see figure (1) as example). 
And in fact experimental analysis of RNA levels in developing embryos shows a wide range 
of gene activity levels. Different genes interact in very complex ways, and activation or 
represion takes place in such a way that the same gene can be more or less active depending 
on the couphng of proteins to different regulatory sites [12]. To see this, let us consider a 
general network formed by m elements showing mutual inhibition, i.e.: 

^ - ...,X.) = - (6) 

With i = 1, ...,m. This network leads to the emergence of a (usually large) number of 
attractors when n is large enough. An example of the phase space trajectories for m = 3 
and II = 1/2 is shown in figure 2. We can see in figure (2) upper that for low n values 
(n = 1) only a coexistence point is reached. However, see figure (2) lower, once the 
nonlinearity is increased beyond a treshold, different attractors are obtained, indicated 
by the flow of the vector field (here n = 4). Three corners are occupied by the simplest 
attractors (molecular switches with P* e {(1, 0, 0), (0, 1, 0), (0, 0, 1)}). But we also have 
P* = (a;*,x*,x*) with x\ ^ x*^ xl e {a,/?, 7} (here a = 0.13, /3 = 0.34,7 = 0.97). So 
three different states are available for each variable. For a more complex network, many 
possible states will be possible displaying a range of intermediate states. In order to avoid 
this problem other approximations have been proposed. These approximations retain the 
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discreteness of the units but consider a continuous range of values [11, 14] more close to 
the real picture, where gene activity can be measured in terms of RNA concentrations 
[12]. 

These results makes necessary a generalization of RBNs extending them into a more 
general framework. In particular, we can ask how arc the critical boundaries defined for 
a more general situation when multiple states are available. 

Although the starting point of our study deals whith genomic regulatory circuits, a set 
of different, but related problems could be explored within this framework. If a given state 
is identified in terms of a cell type, for example, the interactions between diffcrents cells 
(automata) with different cell type will eventually lead to transitions which are known 
to be rather complex and neighbour-dependent [12]. Another example could be a system 
where each unit can have a complex internal state roughly characterized by scalar quantity 
(or "state") which can be, for example, the mean activity. A further scenario where our 
description applies is a system of machines (computers or agents) each one engaged in a 
given task. These machines might be connected in complex ways, and switch from a task 
to another by depending on their actual inputs. 

Discrete dynamical systems involving a given range of states per unit are well known in 
the physics literature. Cellular automata (as described by equation (1)) are just the best 
known example [1]. Another example from statistical mechanics is the q— Potts model 
[13] which generalizes the classical Ising spin systems. It well known that this model 
shows phase transitions with a critical temperature Tc{q) depending on g, the number of 
possible spin states. 

Our goal in this paper is to explore the properties of the critical points for S-state 
random networks. In this context, the analysis of transition phenomena in CA shows 
that, for large r-neighbor hoods, a sharp (second order ) transition occurs for well-defined 
A-values. However this result cannot be traslated to the genome because typically gene 
regulation is a non-local phenomenon [8,12]. We should expect, however, to find some 
agreement between the CA results and those derived from networks involving multiple 
states and random connectivity. 

The paper is organized the follows manner: in section 2, random networks with multi- 
ple states are define and the critical surface in the {K, S, p) space is derived from Derrida's 
annealed approximation. In the the Appendix I associated with this section a maximum 
entropy variational approach is used in order to compute the expected distribution of 
states at criticality. In section 3 the relationship with CA models for large number of 
states S is analysed. In section 4 a variant of Flyvbjerg's method of stable core com- 
putation is introduced in order to obtain an order parameter (the self-overlap) for the 
transition. In the Appendix II associated with this section we compute a Lyapunov ex- 
ponent for the system. Finally, a discussion of the general results and possible extensions 
is given in section 5. 
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2 The model and Derrida's annealed approach 



The previous definition of random Boolean network, where S = {0, 1} can be generahzed 
to a wider set of states, S = {0, l,...,.? — 1} that is to say, to random networks with 
multiple states (RNS). As with RBN, each element cxj (i = 1, A^) receives K inputs; so 
we have: 



where the superindex s of the functions A? indicate that its inputs (and outputs) can take 
S values and the functions A| are randomly chosen from a set S^.- 

As with RBNs, this model exhibits three characteristic dynamical regimes, which we 
will explore. Three examples are shown in figures (3-5) for a A^ = 100 network with 
S = b states, K = 2 and different p-values. Here p is an additional which was introduced 
by B. Derrida et al. [15] in order to make the RBN transition smooth and continuous. 
The parameter p, known as the bias of Aj, is defined as the probability of having as 
output in the function Aj. This parameter is easy to extend in natural form for the RNS 
(see below). The plots correspond to the chaotic (p=0.60), critical (p=0.689) and frozen 
(p=0.8) regimes. 

Now in order to characterize the critical properties of our system, we will applie the 
well known Derrida's method, based on the annealed approximation [15]. This approach 
is a way to avoid the dynamical correlations which appear as the system evolves in time. 
In the thermodynamic limit the original (quenched) model and the annealed counterpart 
share the same phase transition curves [16]. 

We start with two randomly chosen configurations 



which are also randomly taken from the set Cs{N) of all the possible A^-strings (clearly 

'\\Cs{N) = S^). Following Derrida's method, it can be shown that the overlap ai2(t) G 
[0, 1], defined as the normalized number Nauit) of elements with common states in Ci{t) 
and C2{t), will evolve in time following a nonlinear one-dimensional map. In short, let 
Nai2{t + 1) the net overlap after one iteration of Ci{t) and C2{t) under (7). Then a 
new set of connections and A| functions is again taken from Sic, and a new iteration is 
performed over Ci{t + 1) and C2{t + 1). 

For our S'-state random network and assuming that ^ is the probability that two 
outputs of an arbitrary function A| are identical, following the arguments in [15], we have 
a nonlinear map for ai2 given by: 



(7) 



Ci(t)^(a«(t),...,a«(t)) 
C2{t)^{a?it),...,a'^^\t)) 



(8.a) 



(8.6) 




(9) 
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where: a{^(t) is the fraction of elements with identical inputs and [1 — ai|(t)]/'S' is the 
probability for two automata, with at least one different input, being equal ai t + 1. Now 
the critical point where the phase transition takes place (separating the so called frozen 
and chaotic phases) is obtained from the stability condition: 



dai2{t + l) 



dauit) 



a*=l 




< 1 



(10) 



or in other words the following relation between the average network connectivity and the 
number of available states is reached: 



K 



1 

s 



(11) 



We see that for S — 2 (Boolean net) we recover the well known critical point Kc — 2. For 
^ 2 we can approximate the marginal stability relation to: 



K 



1 



;i2) 



We also see that when 5" — > oo the average connectivity moves to 1 (i.e. the critical 

point is reached only at the lowest connectivity). The previous inequality (10) is related 
with the high-temperature value of damage spreading [23] of the Potts model [13]. 

In RNS it is possible to extend Derrida's p-parameter definition: p is defined as the 
probability of having a as output of the function A| (0 is the quiescent state) and the 
other states are equally likely to be present. Although this choice seems too limited, it 
can be justified from variational arguments (see Appendix I). Now the overlap will evolve 
following the ID map: 



au{t + l) = af^{t) + 



^ S-1 



4(t)) 



(13) 



where we have: 



• a{^(i) is the fraction of elements with identical inputs at time t 

• (1 — a^2{t))P^ is the probability of two automata with at least one different input, 
being equal to "0" at t + 1. 

• (1 — a{|(t))(l — pY/{S — 1) gives the probability of finding two automata with at 
least one different input, being equal (but differents from 0) at t + 1. 
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Now, the stability analysis for a^2 = 1 gives us a general condition for the critical 
points: 



Defining a critical surface Pc — Pci^, S) separating the frozen (ordered) phase from the 
chaotic one. 

Many possible characterizations of the two phases and the critical surface can be used. 
In figure (6a-c) we show three examples of the magnetization m (for different RNS with 
N — 100, K — 2 and S — 5 states). Here m is defined as the number of times n a given 
automaton Si is such that Si = over T steps i. e. m = n/T. Then P(m) is a histogram 
giving the distribution of m for the whole system i. e. it counts how many automata 
with the same average magnetization m. We can see, as expected, how P(m) moves from 
a singular distribution with few peaks to a more continuous distribution at the chaotic 
regime. 

A numerical characterization fo these three regimes can be easily obtained by comput- 
ing the frequency distribution N{T) of cycles of different length T. Since these dynamical 
systems are finite, with a maximum number of S^ states, the previous equations (7) lead 
to two types of attractors (cycles and steady states). However, the presence of a criti- 
cal boundary separating the two regimes is made clear by the existence of well-defined 
scaling laws. In figure (6d-f) we show the frequency distribution of periods for the three 
regimes. At the chaotic phase, long-periodic orbits arc rather common. This is due to 
the exponential increase of attractor lengths in this phase. The frozen regime shows the 
opposite tendency: the periods are very short, with an exponential decay with T. At the 
critical boundary, both short and long periods are present, in such a way that in spite of 
the frequent observation of short-period orbits, long-periodic orbits of length T fa 0{N) 
are all ways found. 

In random Boolean networks, the previous observations have been deeply explored 
through both numerical and theoretical approximations [6,7]. A specially important re- 
lation is the dependence of the average period < T > with the system size N. When 
RBN are used with p = 0.5, it is found that the critical phase {K = 2) is characterized 
by a scaling < T >fti \/N and an exponential growth < T{K > 2) >fti 2^^ is obtained 
at the chaotic phase, with B > 1 [2]. We have explored this scaling behavior in the 
RNS counterpart. Two examples of our calculations are shown in figure 8(a-b). The 
average period length has been calculated in two different cases, both varying p (a) and 
the number of states S (b). we can see how the scaling < T{N) >^ N'^ is present at 
criticality and how an exponential increase in < T > also appears as we enter into the 
chaotic phase. Numerical simulations show that for the chaotic phase an exponential law 
< T{S,N) >Ki S^'^ is obtained, with B > 1. Numerical estimations of this value gave 




(14) 



B = 1.11 ±0.02. 
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3 CA and RNS 



It can easily be checked from (14) that for 5" = 2 (Boolean net) we have 




(15) 



which gives us the well known critical line in the (i^T, p) -phase space [15,17]. 

Equation (15) is reminiscent of those found in standard CA (with nearest-neighbor con- 
nectivity [18]). The existence of a transition domain for CA models has been characterized 
by means of several approaches. Langton's A-parameter gives a rough characterization of 
these critical points when complex CA rules (involving high S and/or r values) are used. 
If is the total number of (2r -|- l)-configurations in the CA rule table, as described 

by (1), and M out of the total map to a non-zero state, then A is defined as [18]: 



Thus, 1 — A corresponds to p under the RNS definition. This parameter has been shown to 
be related with -but not equivalent to- a temperature in statistical physics [18]. It was ob- 
served (particularly in the large-r limit) that these CA exhibit some of the characteristics 
of second-order phase transitions (for a detailed study see [18], although the association 
between computation, complexity and critical phenomena has been shown to be fiawed 
[19]). Li et al. [18] showed by means of mean field approximations that, for S = 2, there 
is a simple relationship between the critical value Ac and r (the neighborhood radius): 



This result is consistent with equation (15). But in (17), instead of 2r -|- 1 (the total 
connectivity of the CA model) we have r + 1 (half of it plus one). In the RNS counterpart, 
we have K, the total connectivity. These differences are due to the spatial correlations 
intrinsic to CA models, which are destroyed in the RNS counterpart. This is consistent 
with the analysis of Li et at. [18] which is based in a mean- field estimation of the spreading 
rate of diferent CA patterns. Assuming symmetry in the average spreading rate 7(r) in 
right and left directions, they obtain the onset to non-zero spreading at Ac given by (17). 
The corresponding critical curves for both the 1 — A parameter (squares) and p (circles) as 
a function of the connectivity are shown in figure (8). Remember that a neighbor radius 
r in CA is equivalent to X = 2r -|- 1 in RNS. Thus the first point represented for RNS 
whit 5* = 2 is = 3 with p^. — 0.79. Both curves share the same qualitative behavior, 
but larger values of the bias are required in RNS in relation with CA in order to reach 
the ordered regime, as expected. A similar relationship is found for the temperature in 
relation with q in some Potts models [13]. 




(17) 
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4 Stable core and the self-overlap method 



A well-defined critical phase transition should be characterized by an appropiate order 
parameter Q such that O > at the disordered regime and zero otherwise. In fact, 
d* = 1 — a* as like parameter of order for RBN and RNS [15]. In this context, Flybvjerg 
[22] explored a different approach than Derrida's by defining the stable core at time t, c{t), 
as the (normalized) fraction of automata Si (independent of the initial condition) reach 
stable values at a given step t, i.e. remain constant after t. The stable core asymptotic 
behavior, c{t oo), can be used as an order parameter for the frozen-chaos transition in 
RBN. Explicitly, if f2 = 1 — c(t — > oo), the previous requirements for the order parameter 
definition hold. 

The argument of Flybvjerg for finding an iterated equation for c{t) is as follows. There 
are K + 1 mutually exclusive reasons for a given automata Si to be part of the stable 
core from step t to t + 1. The probabilities of this K + 1 reasons are: (0) All the inputs, 
c"ii(^), C"i2(^); CTjif (^); bcloug to the Stable core at t. Given that the inputs are choosen 
at random, this situation occur with probability c^{t) (where c{t) is interpreted as the 
probability of belonging to the stable core in t) . 

(1) Aj is independent of one input and the rest of the K — 1 variables belong to the 
stable core: 



where pi is defined as the probability of Aj to be independent of one input after fixing its 
K — 1 remainder inputs. 

(j) Aj is independent of j inputs and the rest of the K — j variables belong to the 
stable core: 



where pj is defined as the probability of Aj to be independent of j arguments after fixing 
its K — j remainder arguments. 

{K) A| is a constant fuction (independent of its K arguments), with all out the stable 
core: 




(18) 




(19) 



(1 - c{t)rpK 



(20) 



Adding the K + 1 exclusive reasons: 




(21) 
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where pi is the probabihty that the Boolean function is independent of a certain number 
i of inputs after fixing its X — i remainder inputs (a biologically relevant case are the 
canalizing functions that depend on a unique input and are independent of remainder 
[2]). For the Boolean functions with bias p, i.e., with value with probability p and value 
1 with probability {1 — p) it gives: 

+ (22) 

By substituying (22) into (21) and analyzing its stability the following transition critical 
curve is obtained: i^2p(l — p) = 1 according to [15,23], reached by different methods. 

How to extend this argument to RNS? Apparently, such an extension is straightfor- 
ward. It seems that, again, there are K + 1 exclusive reasons for a given automata to 
move from outside the stable core to be a member of it. Such a simple translation of the 
previous procedure is not possible, and we can show why by means of a simple example. 
Let us consider a RNS with K = 4 and S = 3, i.e.: 

Al^K{ai„ai„ai„ai,) ; i^l,2,...,N (23) 

and let us assume that at a given t the first three inputs belong to the stable core , i.e: 
(Tjj , (Tj2 , CTjg G s{t) and (jj^ ^ s{t). The question is whether or not (Tj(t + 1) e s{t + 1) 
if (Ti{t) ^ s{t). This will happen if Af is independent of the fourth input. But it can 
also occur (and this makes a big difference) if can only reach a subset S' C S, say 
S' = {0, 1}. In such case, we have a reduction in the number of degrees of freedom for 
(Ti4 . In this case it can happen that A? is independent of for cTj^ e E' but not if CTj^ = 2 
(for 5* = 2 a reduction of one degree of freedom simply means that the input belongs to 
the stable core). These possible outcomes forces us to weight the probabilities that the 
input units have their degrees of freedom reduced. More generally, the exclusive reason 
(1) implies that A| is independent of a given input ai^. This occurs with probability 
Pi — p'^ + {1 — pY for the case of two states. For ^'-states, it can occur that Ui. e E' C E 
and Gi. ^ s{t). In order to calculate the order parameter equation for networks with 
multiple states, we will compute the successive overlap (or self-overlap [25]) at+i,t, defined 
as the fraction of automata such that: (Ti{t-\- 1) = crj(t). This is in fact the overlap between 
the system configuration at t with the next configuration at i + l. If the assymptotic value 
of the stable core is 1 the self-overlap at+i^t{t oo) will be one too. This allows us to 
use the asymptotic value of the self-overlap as an appropiate order parameter and as an 
alternative to the standard overlap equation analysed in section 2. 
The iterated map for the self-overlap will be 
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where we desribe the self-overlap between the set at t + 1 and t as a function of the 
self-overlap for {t — l,t). If we interpret at^t_i as the probability for an arbitrary unit to 
remain in the same state at both t — 1 and t, the term a^_j gives the probability that 
all the inputs of a given unit are the same from t — 1 to t. Obviously (1 — o-t^t-i) ^^c 
probability that at least one of the inputs will be different between t and t — 1. In that 
case, there is still a possibility of remaining in the same state at t and at t + 1 just by 
chance. This is given by (p^ + (1 — pY)/{S — 1). 
The stability analysis now gives: 



da 



da- 



t,t-i 



- K 



< 1 



(25) 



which is easily interpreted from a perturbative point of view [23]: / + is the 

probability of no-propagation for the minimal perturbation (one change of state / to m in 
one input). Then, K 1 — [j? + ^^s-i ) mean number of changes generated by the 

perturbation in one time step. This interpretation allows to define a Lyapunov exponent 
for this system (see Appendix II). 

In figures (9-10) we show (continous lines) the asymptotic behavior of both the overlap 
and the self-overlap for different p and K and = 10 from equation (24). As we can 
see, the dynamical equation for the self-overlap (24) is the same as the standard overlap 
equation (13), although it is conceptually different from it. The self-overlap, as it occurs 
with the stable core, can only be applied to the quenched system and fails with the an- 
nealed case. Besides, the self-overlap is computationally more efficient than the standard 
overlap, which requires the parallel running of the two replicas of the system. 

In figure (9) we have followed the time evolution (black symbols) of overlaps between 
two RNS quenched replicas of size N = 1000 as in [15]. Each point is compute averaging 
100 experiments. The connectivity K = 2 and the number of avalible states 5" = 10 are 
fixed, and the bias p change. The white symbols represent identical averages for self- 
overlap, except for the crosses where N = 10.000 for critical p = 0.70 in order to show the 
finite size effect of the critical transition. This effect is more acused for the self-overlap 
as for the overlap. 

In figure (12) we show the time evolution, this time fixint the states 5" = 10 and the 
bias p = 0.81, of the overlap and self-overlap. Again the continous lines are the theoretical 
evolution from equation (29). Black and white symbols correspond to overlap and self- 
overlap respectivily. In this case the variation of the connectivity despite the transition 
ai K — 3. Again the crosses (A^ = 10.000 for self-overlap) show the finite size efecfct. 

All our simulations show a very good agreement whit the theoretical evolution of 
the overlap (eq.l3) and the self-overlap (eq.24). Thus, any of them acts as good order 
parameter. 
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5 Summary 



In this paper we have proposed a simple, straighforward extension of random Boolean 
networks in terms of a discrete model with multiple states. The starting point of our 
analysis is the observation that several natural systems involving a wired set of entities of 
some kind often display a range of state values instead of simple on -off characteristics. In 
some cases, as the specific tasks (states) performed by ants in a colony, the discreteness 
of the states is rather well defined, although their transitions (as a consequence of inter- 
actions of ants engaged in different tasks) can be rather complex [24] suggesting that the 
underlying task-dependent rules are complicated ones. When cell types are considered, 
we also get a well-defined set of " states" although each one itself is the result of a complex 
dynamical system (the gene expression pattern) at a given stationary state. A classical 
set of experiments with different organisms at different levels of development reported 
the existence of complex interactions between cell types eventually leading to switches by 
depending of the cell types under interaction [2,12]. In other systems, like the genomic 
regulatory network, the activity of some genes is sometimes close to an on-off switch but 
it displays a range of activity levels. A first approximation to such continuous networks 
is the RNS described and analysed in this paper. Following previous studies, our interest 
was focused on to the critical features of these networks. It has been conjectured that a 
wide range of complex systems, from evolving ecosystems, biological regulation networks, 
traffic or ant colonies [4] , to cite only a few, display patterns in space and time suggesting 
that they are close to phase transition points. In this context, we have found the explicit 
form of the transition domains as well as an order parameter equation for the transition. 

The critical boundaries for this model have been obtained from Derrida's annealed 
approximation and the computation of s{t) (the stable core equation) required the devel- 
opment of an alternative treatment to the Flyvbjerg derivations for standard RBN [22]. 
This model has been shown to share some common traits with other related dynamic 
models like cellular automata and the g-Potts model. 

Several extensions of this work are possible. Further quantities can be defined in order 
to characterize the damage spreading. This can be easily done by means of theoretical 
extensions of previous definitions of Lyapunov exponents for RBNs [25] . One of them is the 
analysis of the square-lattice counterpart of the RNS with nearest-neighbor connectivity. 
This type of network has been proved to be extremely useful in order to display the critical 
features of the RBN model in terms of percolation on a square lattice [15]. A different 
extension would be a re-definition of the network rules in order to make the transitions 
from different states closer to the continuous dynamics. If the output of the A| function is 
randomly chosen from S, we must expect to observe a time evolution of individual units 
in terms of jumps from a given state to any other state. But if we want to retain the 
relationship between the yS'-state model with the corresponding continuous counterpart 
(as described by sets of nonlinear differential equations) a consistent increase or decrease 
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in individual activity should be observable. This problem can be easily solved through 
the introduction of a new set of rules describing the outputs of the A|-function in terms 
of increases/decreases of gene activity and will be reported elsewhere. 
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6 Appendix I: RNS and maximum entropy 

The relation (11) can be easily generalized: K can be not unique, but an average connec- 
tivity < K > is always able to be defined [17] and we can have a probability distribution 
of states, i.e. P{^) = -P[Af = /i] for ^ E T,. In this general case, the critical line tran- 
sition is given by: < K >= 1/(1 — Y,P'^{l^)) [23]. This result can be used in order to 
find the expected probabihty distribution of states at criticality [17]. Using the following 
maximum entropy constraints: 

Pp*(/x) can be derived from a variational procedure known as the maximum entropy for- 
malism [20]. We first construct the Lagrangian: 

>C(P) = iogP(/^) - p{^P\^A - - «(E^(/^) - i) (3) 

fX IX fX 

where P = (P(0), P(l), . . . , P{S - 1)) and the variation 

is performed. This leads to a set of equations 

logP(//) + 2P(//)/3=-(l + Q;) (5) 
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which must be solved for Let us introduce P{0) = w and P(/i 7^ 0) = a^w with 

Ijl— 1, 2, . . . , 5" — 1 with w e [0, 1] and {ct/i} is a set of positive constants such that 

^1 + E"m) = 1 (6) 

Then we get: 

logw + 2Pw — —(1 + a) 
logw + loga^ + 2Pa^w ^ -{1 + a) = 1, 2, . . . , 5 - 1 (7) 
So = C (constant) V/x and as a consequence 

P(0) = -u; (8.a) 

^(^) = ^ /.= 1,2,...,5-1 (8.6) 

The most hkcly probability distribution close to criticality will be a delta-shaped -P(/^) 
with a given state present with probability w and the rest equally distributed. We should 
note that this property (i.e. a sharply peaked distribution with one or a few predominant 
states) has been reported for the onset of chaos in low-dimensional chaotic dynamical 
systems [21]. Here the probability distributions maximize a previously defined complexity 
measure C and numerical simulations of RNS indeed show that such singular distribution 
with one or a few peaks is characteristic at criticality. 



7 Appendix II: Lyapunov exponents for RNS 

If, in the distance method or Derrida's aproximation, we interpret the replica of the one 
system as a perturbation on the original system [25], we can define a expansion rate of 
the perturbation in a time t for RNS easilly: 



rj{t) 



d{t + l] 



(1) 



In a natural way, we can define now the Lyapunov exponent as the temporal average of 
the logarithm of the expansion rate: 



A(T) = -i:iog^(t)=log 



t=i 



l[v{t) = logr] 



(2) 
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If we use the equation for the distance between configurations: 

di2{t + 1) = (p' + [1 - (1 - d^,{t)f] 

Then the expansion rate in the time t of the perturbation, is: 



(3) 



vit) 



dl2(t+l) 



and approximating for small d(t): 



'2-){i-[i-d(tr} 



1 - d{t)f 1 - Kd{t) 



we have: 



S-1 ^ 

an average constant expansion rate that give us that Lyapunov's exponent: 

{1-pr 



X — lo£ 




S-1 



K 



(4) 

(5) 
(6) 

(7) 



that determines the two classics regimes: A < (order) and A > (chaos), whit the 
marginal case A = 0, in total concordance with the transition surface. 

There are diferent methods to compute the largest Lyapunov exponent. In order to 
show consistence we will demonstrate that it is possible compute Lyapunov exponents 
from the self-overlap in agreement with the previous result. The Wolf's method [26] is 
used to estimate numerically Lyapunov exponents from time scries. In short, the method 
is described as follows: get two points of the temporal serie, let us say X(ti) and X(t2) 
and calculate their distance: | X(t2) — X(ti) |. Assume that | X(t2) — X(ii) |< e, being 
e > small. Next, compute the distance, after a time T, i.e: | X(t2 +T) — X(ti + T) |. 
This time, usuallyly, is a fraction of the characteristic period or the time required for the 
autocorrelation function go to zero. Repeating for n pairs of points and averaging, we 
obtain an estimation of the Lyapunov exponent: 



A = 4; E log 



1 X(t2 + T) 


-X(ti+T) 1 


1 X(t2) 


- X(ti) 1 



(8) 



For RNS, we can write down an equation for the normalized Hamming distance between 
succesive time steps in our system, i.e. for the complementary probability of the auto- 
overlap: dt^t-i — ^ — cit,t-i- Thus the equation for the self-overlap becomes: 



d 



t+i,t 



1 - 




S-l 



[1 - (1 - 



(9) 
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If we approximate linearly close to the fixed point d* = 0. The function becomes: 



dt+i,t = K 
The iterated equation now gives: 



S-l 



t,t-i 



(10) 



t+T,t+T-l 



K 



(11) 



To compute (8) in our terms, and taking t2 — ti-\- 1, i.e.: a natural step in our system, 
we have: 



d-, 



t+T,t+T-l 



dt,t-i 



K 



S-l 



(12) 



a constant vahic that, once introduced in the sum (8), determines the Lyapunov exponent 
(7) as sugcstcd previously. 
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Figure 1: Gene network with N = 3. Here the proteins (PI, P2, P3) obtained from each 
gene interact in complex ways with DNA. The gene products bind some specific sites 
either in the DNA sequence or on another protein site, thus inhibiting or stimulating gene 
activity. Several genes can be involved in the regulation of a given gene (as described 
in RBN models through the connectivity K. But such regulation can have different 
effects on different genes and under different initial conditions. Here genes 2 and 3 act 
synergistically stimulating gene 3 in such a way that the latter produces at a rate twice 
the one expected if gene 2 were repressed by gene 1. If gene 1 starts producing an enough 
large concentration of PI, then it would be able to supress the stimulation from gene 2 to 
3 and thus gene 3 would produce a basal level of protein. If we identify such differences of 
gene activity levels as gene states then wc can map the previous example into a discrete 
model with a number of (multiple) states 



Figure 2: Three-dimensional vector fields for the nonlinear system (6) for m = 3, /x = 1/2 
and two different values of the nonlinearity n: n — 1 (upper diagram) and n — A (lower). 
For a small degree of nonlinearity a single attractor is present, but an increase in the 
strenght of the interactions leads to a multiplicity of attractors, both in some of the 
corners of the cube and inside it 



Figure 3: Space-time diagram for a chaotic RNS, with N — 100, S = 5, K — 2 and bias 
p = 0.60 
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Figure 4: Space-time diagram for a RNS at the critical boundary, with N — 100, S — 
5, X = 2 and bias p = 0.689 

Figure 5: Space-time diagram for a RNS at the frozen regime, with N — 100, S — 5, K — 2 
and bias p = 0.80 

Figure 6: (a-c) distribution of magnetization m for the three different regimes of a RNS 
with K = 2, S = 5. Here: (a) p = 0.60, chaos, (b) p = 0.689 , critical and (c) p = 0.8, 
frozen regime, respectively. These distributions have been calculated using a, N — 100 
network and averaging over r — 10^ steps after lO'* transients have been removed. In 
(d-f) the corresponding distributions of cycle lengths are shown for the same parameter 
combinations. M = 10^ different RNS have been used, with random connectivities and 
random initial conditions 

Figure 7: Average period length for RNS with (a) S — 5, K — 2 and different biases and 
for (b) p — 0.5, K — 2 and different numbers of states 

Figure 8: Phase diagram obtained from mean field theories for both RNS and CA models 
{S = 2). Here the fraction of transitions towards the quiescent state is plotted against 
neighbor radius 

Figure 9: Continous lines: dynamical evolution of the self-overlap of RNS of K = 2 and 
5* = 10 with initial self-overlap (and overlap) of 0.1 and differents values of the bias p 
through the iteration of equation (24). The black symbols shows numerical averages of 
the overlap between two RNS quenched replicas for 100 experiments with RNS of size 
N — 1.000. The white symbols shows identical averages and conditions for self-overlap, 
except for the crosses where = 10.000 for critical p — 0.70 in order to show the finite 
size effect of the critical transition 

Figure 10: Continous lines: dynamical evolution of the self-overlap of RNS of p = 0.81 
and 5" = 10 with initial self-overlap of 0.5 and differents values of the connectivity K 
through the iteration of equation (24). The black symbols shows numerical averages of 
the overlap between two RNS quenched replicas for 100 experiments with RNS of size 
A^ = 1.000. The white symbols shows identical averages for self-overlap, except for the 
crosses where A^ = 10.000 for critical connectivity A' = 3 in order to despite the finite 
size effect of the critical transition 



20 



0.2 0.4 0.6 0.8 
Magnetization 



0.2 0.4 0.6 0.8 1 




30 



o 20 

CD 

.? 10 











3000 



o 2000 

CD 



CT 
CD 



1000 











8000 



>, 6000 

o 
c 

^ 4000 

0) 

^ 2000 











p=0.60 
Chaotic 




1000 



2000 



p=0.68 
Critical 



50 



p=0.8 
Frozen 



10 
Period 



100 



20 




0.6 ^ ' ^ ' ^ ' ^ ' ^ ' ' 

10 20 30 40 50 

t 




0.0 ' ' ^ ' ^ ' ^ ' ^ ' ' 

10 20 30 40 50 

t 



